function showFamilyOfTraces(varargin)
% showFamilyOfTraces
%
% A userFcn for qvTrace
%
% Runs from any trace within the family.

% gs 2008 03 08
% gs 2009 05 30 -- make it work with data acq'd with older ephus versions
% ben 2009-09-03 -- prompt for holding condition (resting, etc) and prefix
% the output structure appropriately
% 2009-10-07 Ben Suter -- handle case where not enough spikes to compute
% some fI parameters: set to NaN to avoid error
% 2009-10-07 Ben Suter -- add VsagpercentAlt, an alternative sag
% calculation (see below for more detail)
% 2009-10-07 Ben Suter -- plot V-I curve (figure has 5 rows now and is
% taller), and calculate input resistance as slope of a linear fit to the
% V-I curve.
% 2009-10-08 Ben Suter -- add VsagAltPeakTimeMs, the time to peak in ms
% since step onset. TODO: A more important measure is when the response
% reaches steady-state, but I have not yet implemented this. Also renamed
% VsagpercentAlt to VsagAltPercent, for consistency.
% 2009-10-12 Ben Suter -- added V_baseline calculation, mainly for resting
% potential, also to confirm ca. -70 mV baseline when holding; also removed
% the "note" line since I don't use it.
% 2010-02-02 Taro Kiritani -- added channel num mode for multipele channels.
% ---------------------------------

% LOAD TRACE DATA
currentTrace = evalin('base', 'currentTrace');
sr = currentTrace.header.ephys.ephys.sampleRate;
Ntot = currentTrace.header.loopGui.loopGui.iterations;
Nseq = currentTrace.header.loopGui.loopGui.iterationCounter;
Nacq = str2double(currentTrace.header.xsg.xsg.acquisitionNumber);
[p, f, ext] = fileparts(currentTrace.userFcn.filename);
if strcmp(currentTrace.userFcn.fIAnalysis.familyType, 'cpu loop')
    NfirstFile = Nacq - Nseq;
    % happens to come out right; Nacq actually seems to be for
    % the next, i.e., to-be-acq'd, trace
    for n = 1 : Ntot
        newNum = num2str(NfirstFile + n - 1);
        nchar = numel(newNum);
        f(end-3:end) = '0000';
        ftemp = [f(1 : end - nchar) newNum];
        filename = fullfile(p, [ftemp, ext]);
        currentTrace = load(filename, '-mat');
        
        if nargin > 0
            if varargin{1} == 1
                traceData(:,n) = currentTrace.data.ephys.trace_1;
            elseif varargin{1} == 2
                traceData(:,n) = currentTrace.data.ephys.trace_2;
            else
                traceData(:,n) = currentTrace.data.ephys.trace_1;
            end
        else
            traceData(:,n) = currentTrace.data.ephys.trace_1;
        end
        
        traceLength = currentTrace.header.ephys.ephys.traceLength;
        cp = currentTrace.header.pulseJacker.pulseJacker.currentPosition;
        % conditional; for backward compatibility; gs 2009 05 30
        if(currentTrace.header.xsg.xsg.version) >= 0.2
            if nargin > 0
                if varargin{1} == 1
                    params = currentTrace.header.pulseJacker.pulseJacker.pulseDataMap{1, cp + 1};
                elseif varargin{1} == 2
                    params = currentTrace.header.pulseJacker.pulseJacker.pulseDataMap{2, cp + 1};
                end
            else
                params = currentTrace.header.pulseJacker.pulseJacker.pulseDataMap{1, cp + 1};
            end
        else % older traces
            np = size(currentTrace.header.pulseJacker.pulseJacker.positions{1}, 2);
            params = currentTrace.header.pulseJacker.pulseJacker.pulseDataMap{(cp-1)*np+1}(:);
        end
        pulseTrace = pulseTraceFromParameters(params, traceLength);
        Isteps(n) = params.amplitude;
        %     timedata = ((0 : currentTrace.header.ephys.ephys.sampleRate * currentTrace.header.ephys.ephys.traceLength - 1) / ...
        %         currentTrace.header.ephys.ephys.sampleRate)';
        stimData(:,n) = pulseTrace;
    end
elseif strcmp(currentTrace.userFcn.fIAnalysis.familyType, 'daq loop')
    filename = fullfile(p, [f, ext]);
    traceLength = currentTrace.header.ephys.ephys.traceLength; % in seconds
    for n = 1 : Ntot
        startInd = 1 + (n-1) * sr * traceLength;
        endInd = n * sr * traceLength;
        if nargin > 0
            if varargin{1} == 1        
                traceData(:,n) = currentTrace.data.ephys.trace_1(startInd:endInd);
            elseif varargin{1} == 2
                traceData(:,n) = currentTrace.data.ephys.trace_2(startInd:endInd);
            else
                traceData(:,n) = currentTrace.data.ephys.trace_1(startInd:endInd);        
            end
        else
            traceData(:,n) = currentTrace.data.ephys.trace_1(startInd:endInd);        
        end
%         cp = currentTrace.header.pulseJacker.pulseJacker.currentPosition
        np = size(currentTrace.header.pulseJacker.pulseJacker.positions{1}, 2);
%         params = currentTrace.header.pulseJacker.pulseJacker.pulseDataMap{(cp-1)*np+1}(:);
        params = currentTrace.header.pulseJacker.pulseJacker.pulseDataMap{(n-1)*np+1}(:);
        pulseTrace = pulseTraceFromParameters(params, traceLength);
        Isteps(n) = params.amplitude;
        %     timedata = ((0 : currentTrace.header.ephys.ephys.sampleRate * currentTrace.header.ephys.ephys.traceLength - 1) / ...
        %         currentTrace.header.ephys.ephys.sampleRate)';
        stimData(:,n) = pulseTrace;
    end
else
    disp('showFamilyOfTraces: currently only works with cpu or daq loops');
    return
end


% FIGURE
figure('Name', [f(1:6) ' family'], 'NumberTitle', 'off', 'Color', 'w', ...
    'Position', [264    36   940   900]);
colormap(flipud(gray(256)));
xdata = (0:size(traceData, 1)-1) / sr;
Nrows = 5; Ncols = 3; plotnum = 0;

% CURRENT STEPS
plotnum = plotnum + 1;
subplot(Nrows, Ncols, plotnum);
if max(size(xdata)) ~= max(size(stimData))
    stimData = stimData(1:max(size(xdata))); % trim
end
plot(xdata, stimData); 
box off; ylabel('Current (pA)');

% TRACES, OFFSET
[samples, numTraces] = size(traceData);
traceData2 = traceData + repmat((0:numTraces-1)*50, samples, 1);
plotnum = plotnum + 1;
subplot(Nrows, Ncols, plotnum);
plot(xdata, traceData2);
axis off; axis tight; 

% TRACES, SUPERIMPOSED
plotnum = plotnum + 1;
subplot(Nrows, Ncols, plotnum);
plot(xdata, traceData);
box off; ylabel('Vm (mV)');

% TRACES, AS AN IMAGE
plotnum = plotnum + 1;
imgplot = subplot(Nrows, Ncols, plotnum);
himg = imagesc(traceData');
xlims = [xdata(1) xdata(end)];
set(himg, 'XData', xlims, 'YData', Isteps);
set(gca, 'CLim', [-100 0], 'YDir', 'normal');
axis tight; box off; xlabel('Time (sec)'); ylabel('I step (pA)');

% SPIKE DETECTION
delay = params.squarePulseTrainDelay;
dur = params.squarePulseTrainWidth;
for N = 1 : numTraces
    stepData = traceData(delay*sr : (delay+dur)*sr, N)';
    traceFilt = stepData - medfilt1(stepData, 100);
    a = traceFilt > 5; % NB: thresholding set here
    spikeRight = find(diff(a) < 0);
    spikeLeft = find(diff(a) > 0);
    n = length(spikeRight);
    numberOfSpikesPerStep(N) = n;
    if n == 0 
        spikeIndex = NaN;
        spikeAmplitude = NaN;
        spikeTimes = NaN;
        meanSpikeRate(N) = n;
    else
        spikeIndex = []; spikeAmplitude = [];
        for i = 1 : n
            [Y I] = max(stepData( spikeLeft(i) : spikeRight(i) ));
            spikeIndex(i) = I + spikeLeft(i) - 1;
            spikeAmplitude(i) = stepData( spikeIndex(i) );
            spikeTimes = spikeIndex / sr;
            meanSpikeRate(N) = round(n / dur);
        end
    end
    if n >= 2
        interSpikeInterval = ( spikeTimes(2) - spikeTimes(1) );
        instRate1and2(N) = round(1 / interSpikeInterval);
    else
        instRate1and2(N) = NaN;
    end
    spikedata(N).spikeTimes = spikeTimes;
    spikedata(N).spikeAmplitude = spikeAmplitude;
end
axes(imgplot); hold on
for n = 1 : numel(spikedata)
    plot(spikedata(n).spikeTimes+delay, Isteps(n), 'ro');
end

% assignin('base', 'spikedata', spikedata);


% TODO: if no spikes at all, insert null values and exit

% MEAN FREQUENCY VERSUS CURRENT
plotnum = plotnum + 1;
subplot(Nrows, Ncols, plotnum);
plot(Isteps, meanSpikeRate, 'bo-'); 
hold on 
xlabel('I (pA)'); ylabel('f (Hz)');
if any(meanSpikeRate) > 0
    inds = find(meanSpikeRate > 0); % find traces with spikes
    Ithresh1Xind = inds(1); % find the threshold I step
    Ithresh = Isteps(Ithresh1Xind);
    Ithresh2Xind = find(Isteps==Ithresh*2); % find the 2X I step (may not exist)
    if ~isempty(Ithresh2Xind)
        indnew = find(inds==Ithresh2Xind);
        inds = inds(1:indnew); % exclude steps above 2X threshold
    else
        Ithresh2Xind = numel(Isteps); % defaults to the last I step
        disp(' '); disp('2X-threshold step not found; using all available steps.')
    end
    inds = [inds(1)-1 inds]; % include the last zero-spike Istep as well
    p = polyfit(Isteps(inds), meanSpikeRate(inds), 1);
    fIslope = p(1);
    fIoffset = p(2);
    yfit = p(1) * Isteps(inds) + p(2);
    plot(Isteps(inds), yfit, 'm-');
    title(['f-I: mean slope = ' num2str(p(1)) ' Hz/pA']);
    hline = line([Ithresh Ithresh], get(gca, 'YLim')); % line at threshold
    set(hline, 'Color', 'y', 'LineStyle', '-');
    if ~isempty(Ithresh2Xind) % line at 2X threshold
        hline = line([Ithresh*2 Ithresh*2], get(gca, 'YLim'));
        set(hline, 'Color', 'y', 'LineStyle', '--');
    end
else
    Ithresh2Xind = [];
    title('f-I (no spikes)'); fIslope = NaN; Ithresh = NaN;
end
if any(instRate1and2) > 0
    inds2 = intersect(inds, find(instRate1and2 > 0)); % traces with >=2 spikes
    p2 = polyfit(Isteps(inds2), instRate1and2(inds2), 1);
    fIslopeFirstISI = p2(1);
    fIoffsetFirstISI = p2(2);
    yfit2 = p2(1) * Isteps(inds2) + p2(2);
    plot(Isteps, instRate1and2, 'co-');
    plot(Isteps(inds2), yfit2, 'm-');
else
    % Set these to NaN if unable to compute them (Ben Suter)
    fIslopeFirstISI = NaN; 
    fIoffsetFirstISI = NaN;
end

% TRACE WITH MANY SPIKES -- SPIKE FREQUENCY ADAPTATION
% SFA: 3rd/5th ISI, in 1st trace with >=6 spikes
plotnum = plotnum + 1;
subplot(Nrows, Ncols, plotnum);
a = find(numberOfSpikesPerStep>=6);
if ~isempty(a)
    traceInd = a(1);
else
    traceInd = numel(Isteps);
end
plot(xdata, traceData(:, traceInd));
latencyvec = spikedata(traceInd).spikeTimes;
intervalvec = diff(latencyvec);

if length(intervalvec) >= 5
    SFAratio = intervalvec(3)/intervalvec(5);
else
    disp('Unable to compute the ratio');
    SFAratio = NaN;
end
box off; ylabel('Vm (mV)');
if ~isempty(Ithresh2Xind)
    r = Isteps(Ithresh2Xind) / Isteps(Ithresh1Xind);
    title(['SFA: 3rd/5th ISI = ' num2str(SFAratio)]);
end
% SFAratioIstepAmp = r; % -fold amplitude of Istep over threshold

% INTERSPIKE INTERVALS -- SFA, cont.
% Convert spike times to interspike intervals
for n = 1 : numel(spikedata)
    istep = Isteps(n);
    t = spikedata(n).spikeTimes;
    if numel(t) > 2
        dt = diff(t);
        f = dt.^-1;
        sfa = f(2:end)./f(1:end-1);
    else
        f = NaN;
        sfa = NaN;
    end
    f_mean(n) = mean(f);
    sfa_mean(n) = nanmean(sfa);
end
assignin('base', 'f_mean', f_mean)
assignin('base', 'sfa_mean', sfa_mean)
% INSTANTANEOUS FIRING RATES
% plotnum = plotnum + 1;
% subplot(Nrows, Ncols, plotnum);


% TRACE AT 1X THRESHOLD
plotnum = plotnum + 1;
subplot(Nrows, Ncols, plotnum);
traceAtThresh = traceData(:, Ithresh1Xind);
plot(xdata, traceAtThresh);
box off; ylabel('Vm (mV)'); title('Response at threshold');

% SUPERIMPOSED SPIKES
plotnum = plotnum + 1;
subplot(Nrows, Ncols, plotnum);
% get snippets around peaks
tpre = 0.005; % pre-peak dur, sec
ind1 = round(tpre * sr);
tpost = 0.02; % post-peak dur, sec
ind2 = round(tpost * sr);
tpostLong = 0.12;
indLong = round(tpostLong * sr);
snippetThreshSpike = [];
for n = 1 : numel(spikedata)
    tpeaks = spikedata(n).spikeTimes+delay;
    if ~isnan(tpeaks)
        snippets = [];
        for i = 1 : numel(tpeaks)
            indpeak = round(tpeaks(i) * sr);
            snippets(:,i) = traceData(indpeak - ind1 : indpeak + ind2, n);
            if isempty(snippetThreshSpike) % grab the 1st one for later use
                snippetThreshSpike = snippets(:,i);
                snippetThreshLong = traceData(indpeak - ind1 : indpeak + indLong, n);
            end
        end
        plot(snippets);
        hold on
    end
end
axis off; axis tight;

% AP THRESHOLD, HEIGHT, WIDTH
plotnum = plotnum + 1;
subplot(Nrows, Ncols, plotnum);
plot(snippetThreshSpike, 'b'); % first spike at I thresh
hold on; 
dvdt = diff(snippetThreshSpike); % dV/dT for same
plot(dvdt, 'c-');
axis off; axis tight; title('AP properties at threshold')
xlabel('Samples')
ylabel('mV or dV/dt')
% Find where dV/dt first takes off -- e.g. 10% of max
dvdt_baseline = mean(dvdt(1:10)); % grab some values at start of trace
dvdt_max = max(dvdt); % find peak
dvdt_threshlevel = 0.1; % set to 10%; arbitrary ...
dvdt_threshvalue = (dvdt_max - dvdt_baseline) * dvdt_threshlevel;
dvdt_threshindex = find(dvdt > (dvdt_baseline + dvdt_threshvalue));
% interpolate to estimate threshold-crossing time in dvdt trace precisely
y1 = dvdt(dvdt_threshindex(1)-1);
y2 = dvdt(dvdt_threshindex(1));
% x1 = 0; x2 = 1;
slope = y2 - y1;
y3 = dvdt_threshvalue - y1;
x3 = y3 / slope; % i.e., add this fraction of a time increment
% check that it's a correct fraction 
if x3 >= 0 && x3 <= 1
%     disp(['interpolation OK: ' num2str(x3)]);
else
    disp(['interpolation not OK: ' num2str(x3)]);
end
% now use this to interpolate the voltage trace
extraVm = x3 * (snippetThreshSpike(dvdt_threshindex(1)) - ...
    snippetThreshSpike(dvdt_threshindex(1)-1));
% snippetThreshSpike(dvdt_threshindex(1)-1)
Vthresh = snippetThreshSpike(dvdt_threshindex(1)-1) + extraVm;
VthreshTime = dvdt_threshindex(1) + x3;
% snippetThreshSpike(dvdt_threshindex(1))
plot(VthreshTime, dvdt(dvdt_threshindex(1)) + y3, 'ro');
plot(VthreshTime, Vthresh, 'ro');
% AP HEIGHT: relative to threshold
Vmax = spikedata(Ithresh1Xind).spikeAmplitude(1);
plot(ind1, Vmax, 'ro')
Vheight = Vmax - Vthresh;
% AP HALF-WIDTH: taken at half height (with interpolation)
Vhalfheight = Vthresh + Vheight/2;
Vrange = find(snippetThreshSpike(1 : round(2 * tpre * sr)) > Vhalfheight);
% first, interp rising edge
v1 = snippetThreshSpike(Vrange(1) - 1);
v2 = snippetThreshSpike(Vrange(1));
extrafrac1 = (Vhalfheight - v1) / (v2 - v1);
if extrafrac1 >= 0 && extrafrac1 <= 1
%     disp(['interpolation2 OK: ' num2str(extrafrac1)]);
else
    disp(['interpolation2 not OK: ' num2str(extrafrac1)]);
end
% now interp falling edge
v3 = snippetThreshSpike(Vrange(end));
v4 = snippetThreshSpike(Vrange(end)+1);
extrafrac2 = (Vhalfheight - v4) / (v3 - v4);
if extrafrac2 >= 0 && extrafrac2 <= 1
%     disp(['interpolation3 OK: ' num2str(extrafrac2)]);
else
    disp(['interpolation3 not OK: ' num2str(extrafrac2)]);
end
numSamples = numel(Vrange) + (1 - extrafrac1) + (1 - extrafrac2);
% finally calc the halfwidth
APhalfwidth = (numSamples / sr) / 2;
hline = line([Vrange(1)+extrafrac1 Vrange(end)+1-extrafrac2], ...
    [Vhalfheight Vhalfheight]);
set(hline, 'Color', 'k');

% fAHPs,ADPs: post-spike afterpotentials (at threshold)
% plot for Vm
plotnum = plotnum + 1;
hax1 = subplot(Nrows, Ncols, plotnum);
plot(snippetThreshLong);
box off; axis tight; xlabel('Samples'); ylabel('Vm (mV)'); 
title('Afterpotentials');
hold on
plot(VthreshTime, Vthresh, 'ro')
set(gca, 'YLim', [min(snippetThreshLong)-2 Vthresh+2])

% helper plot for dV/dt calc's
plotnum = plotnum + 1;
hax2 = subplot(Nrows, Ncols, plotnum);
dvdt1 = diff(snippetThreshLong);
% plot(dvdt1, 'k');
w = 13;
dvdt2 = gfilter(dvdt1, w, 'avg');
% plot(dvdt2, 'b');
dvdt3 = gfilter(dvdt2, 2*w, 'avg');
% plot(dvdt3, 'r');
dvdt4 = gfilter(dvdt3, 3*w, 'avg');
plot(dvdt4, 'k');
hold on
box off; axis tight; xlabel('Samples'); ylabel('dV/dt');
set(gca, 'YLim', [-.1 .1]); title('Afterpotential calculations')
xlims = get(gca, 'XLim');
line(xlims, [0 0])
% look for upward zero crossing just after AP peak: pre-ADP trough
overzero = dvdt4 > 0;
% plot(overzero*.05, 'm-') 
plot(diff(overzero), 'y-');
inds0Xup = find(diff(overzero) > 0);
inds0Xdown = find(diff(overzero) < 0);
% TODO: ensure that first upwards crossing is after spike
if numel(inds0Xup) >= 1
    Vfahp = snippetThreshLong(inds0Xup(1));
    fAHP =  Vfahp - Vthresh;
    plot(inds0Xup(1), 0, 'ro', 'Parent', hax2)
    plot(inds0Xup(1), Vfahp, 'ro', 'Parent', hax1)
    if numel(inds0Xdown) >= 2
        Vadp = snippetThreshLong(inds0Xdown(2));
        ADP = Vadp - Vfahp;
        plot(inds0Xdown(2), 0, 'ro', 'Parent', hax2)
        plot(inds0Xdown(2), Vadp, 'ro', 'Parent', hax1)
        if numel(inds0Xup) >= 2
            Vmahp = snippetThreshLong(inds0Xup(2));
            mAHP =  Vmahp - Vthresh;
            plot(inds0Xup(2), 0, 'ro', 'Parent', hax2)
            plot(inds0Xup(2), Vmahp, 'ro', 'Parent', hax1)
        else % must have been a continuing AHP after the ADP
            mAHP = NaN;
        end
    else % must have been an AHP but no ADP
        ADP = NaN; mAHP = NaN;
    end
else % didn't reach end of AHP ...
    fAHP = NaN; ADP = NaN; mAHP = NaN;
end

% "Find the min within 50 ms of the spike, and the max within
% 70 ms of this min. If min is more hyperpolarized than the both the AP 
% thresh and the max, the fAHP is defined as the difference between the 
% AP thresh and the mean Vm in a small time window about the min."
%
% "The DAP is the difference between the fAHP min and the mean Vm in 
% a small time window about the max."

% SAG POTENTIALS
% e.g. Stafstrom et al. (1984)
% sag defined here as Vpeak/Vss (both rel. to Vbase)
% or, for over- or undershoot, as Vtail/Vbase (both rel. to Vss) 
%
% Sag and overshoot for hyperpolarizing step of -100 pA
% calc baseline Vm's over first 100 msec (nearly)
% assumes step delay = 100 ms, dur =500 ms
% search for min in 1st 200 ms after step
% determine Vss over last 100 ms (nearly) of step
% search for max in 1st 200 ms after step termination
del = 0.1; dur = 0.5;
i_neg = find(Isteps < 0);
if ~isempty(i_neg)
    for n = 1 : numel(i_neg)
        Vb = mean(traceData(1:round(del*sr)-20, i_neg(n)));
        Vmin = min(traceData(round(del*sr):round((del+0.2)*sr), i_neg(n)));
        Vminamp = abs(Vmin - Vb);
        Vss = mean(traceData(round((del+dur-0.1)*sr):round((del+dur)*sr)-20, i_neg(n)));
        Vssamp = abs(Vss - Vb);
        Vsagpercent(n) = 100 * (Vminamp - Vssamp)/ Vssamp;
        Vomax = max(traceData(round((del+dur)*sr):round((del+dur+0.2)*sr), i_neg(n)));
        Vomaxamp = abs(Vomax - Vss);
        Vovershootpercent(n) = 100 * (Vomaxamp - Vssamp) / Vssamp;

        % 2009-10-07 Ben Suter -- alternative sag calculation, similar to Shin et
        % al. 2008 Neurobiology of Disease. Finds min value during first
        % 200 ms of step, averages 10 ms around this minimum to get vmax. Compares
        % this to the average response over the last 50 ms of the step.
        % Note: Only uses 1st step in train
        if n == 1
            sta = round(del*sr);
            sto = round((del+dur)*sr)-1;
            [m, mi] = min(traceData(sta:round(sta+0.2*sr), i_neg(n)));
            VsagAltPeakTimeMs = 1000 * mi / sr; % time of peak in ms since step onset
            mi = sta - 1 + mi;            
            vmax = mean(traceData(mi-50:mi+50), i_neg(n)) - Vb; % +/- 5ms of minimum
            vss = mean(traceData(sto-0.05*sr:sto, i_neg(n))) - Vb; 
            VsagAltPercent = 100*(vmax-vss)/vmax;
            
            % TODO: Calculate time at which steady-state is reached
%             snip = traceData(mi:sto, i_neg(n));
%             ws = 100; smoothed = filter(ones(ws,1), ws, snip); 
%             figure; hold on; plot(snip, 'b'); 
%             plot(ws:length(smoothed), smoothed(ws:end), 'r');
        end
    end
end

Vsagpercent = mean(Vsagpercent);
Vovershootpercent = mean(Vovershootpercent);

plotnum = plotnum + 1;
subplot(Nrows, Ncols, plotnum);
hold on
% calc baseline Vm's for each step trace
lag = -0.05; % pre-step time (sec)
ind = round( (delay+lag)*sr );
% plot(Isteps, traceData(ind, :), 'ko-');
Vbaseline = mean(traceData(ind, :)); 
% plot Vm at 100 and 300 msec
lag = 0.1; % when to look, post-step (sec)
ind = round( (delay+dur+lag)*sr );
Vsag_100ms = traceData(ind, :);
plot(Isteps, Vsag_100ms, 'ro-');
lag = 0.3; % when to look, post-step (sec)
ind = round( (delay+dur+lag)*sr );
Vsag_300ms = traceData(ind, :);
plot(Isteps, Vsag_300ms, 'go-');
xlims = get(gca, 'XLim');
ylims = get(gca, 'YLim');
hline(1) = line([0 0], [ylims(1) ylims(2)]);
hline(2) = line([xlims(1) xlims(2)], [Vbaseline Vbaseline]);
set(hline, 'Color', 'k');
xlabel('I step (pA)')
ylabel('Vm');
title('Sag (100 & 300 ms post-step)');

%%%%%% Rn: Input resistance from linear range of V-I %%%%%%%%%%%%%%%%%%%%
% 2009-10-07 Ben Suter -- adding this per Shin et al. (2008)
% Omit traces that have spikes, because difficult to compute the baseline.
% In Shin et al. I believe they use much smaller current injection, such
% that the response is within a +/- 3 mV envelope. Thus there are no spikes
% and the effect of active conductances are minimized.
% I just use all non-spiking responses and fit a straight line to the
% resulting I-V curve.
%
del = 0.1; dur = 0.5;
sta = round(del*sr);
sto = round((del+dur)*sr)-1;
vresp = NaN(numel(Isteps), 1);
for n = 1 : numel(Isteps)
    if numberOfSpikesPerStep(n) > 0 
        % Skip any response with spikes, leave vresp as NaN
        continue;
    end
    % Baseline is all before step
    vbl = mean(traceData(1:sta - 1, n)); 
    % Steadystate is last 50 ms of step
    vss = mean(traceData(round(sto - (0.05)*sr):sto, n));
    vresp(n) = vss - vbl;
end

% Assume V(I) = Rn*I + d and fit accordingly (after removing NaN points)
X = [Isteps' ones(size(Isteps'))];
notNaN = find(~isnan(vresp));
fitCoeffs = X(notNaN, :) \ vresp(notNaN); 
inputResistanceMOhm = 1000*fitCoeffs(1);   % fitCoeffs(1) = slope = Rn

% V-I Curve
plotnum = plotnum + 1;
subplot(Nrows, Ncols, plotnum);
plot(Isteps(notNaN), vresp(notNaN)); 
hold on; 
plot(Isteps(notNaN), fitCoeffs(1)*Isteps(notNaN) + fitCoeffs(2), 'r--');
xlabel('Injected current (pA)');
ylabel('Vss (mV)'); % Steady-state voltage
title('V-I, slope = Rin');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% V_baseline: Baseline membrane potential %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For I-steps without holding current, this is resting membrane potential
%
flatStep = (Isteps == 0);
if flatStep < 1
    error('Unable to calculate V_baseline since not step with 0 pA');
else
    V_baseline = mean(traceData(:,flatStep));
end

% store relevant data in a single structure variable
traceFamily.filenameOfLastTrace = filename;
traceFamily.traceData = traceData;
traceFamily.xdata = xdata;
traceFamily.stimData = stimData;
traceFamily.Isteps = Isteps;
traceFamily.fIslope = fIslope;
traceFamily.fIoffset = fIoffset;
traceFamily.fIslopeFirstISI = fIslopeFirstISI;
traceFamily.fIoffsetFirstISI = fIoffsetFirstISI;
traceFamily.Ithresh = Ithresh;
traceFamily.spikeTimes = spikedata;
traceFamily.SFAratio = SFAratio;
traceFamily.sfa_mean_vec = sfa_mean;
traceFamily.sfa_overallmean = nanmean(sfa_mean);
traceFamily.f_mean_vec = f_mean;
traceFamily.Vthresh = Vthresh;
traceFamily.Vheight = Vheight;
traceFamily.APhalfwidth = APhalfwidth;
traceFamily.fAHP = fAHP;
traceFamily.ADP = ADP;
traceFamily.mAHP = mAHP;
traceFamily.Vbaseline = Vbaseline;
traceFamily.Vsag_100ms = Vsag_100ms;
traceFamily.Vsag_300ms = Vsag_300ms;
traceFamily.Vsagpercent = Vsagpercent;
traceFamily.VsagAltPercent = VsagAltPercent;
traceFamily.VsagAltPeakTimeMs = VsagAltPeakTimeMs;
traceFamily.Vovershootpercent = Vovershootpercent;
traceFamily.InputResistanceMOhm = inputResistanceMOhm;
traceFamily.V_baseline = V_baseline;
assignin('base', 'traceFamily', traceFamily);

% %%%%%%%%%%%% output the results in .m-file syntax %%%%%%%%%%
% %% 2009-09-03, Ben Suter: prompts for prefix to distinguish
% %% current steps at various holding currents (e.g. resting, -70 mV)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine the prefix (name of struct), and optionally the holding current
holdingMode = {'resting', 'minus70mV'};
[Selection,ok] = listdlg('ListString', holdingMode, 'SelectionMode', 'single', 'Name', 'IC: No holding (resting) or at -70mV?', 'ListSize', [250 50]);
if ~ok
    prefix = 'map.isteps.default.';
    holdingCurrent = NaN;
else
    if Selection == 1
        prefix = 'map.isteps.resting.';
        holdingCurrent = 0;
    else
        prefix = 'map.isteps.minus70mV.';
        answer = inputdlg({'Holding current (in pA)'}, 'What was the holding current?');
        if isempty(answer)
            holdingCurrent = NaN;
        else
            holdingCurrent = str2double(answer{1});
        end
    end
end

disp(' ')
disp('% ********** I step family **********')
disp([prefix 'holdingCurrent = ' num2str(holdingCurrent) ';']);

numSteps = size(traceFamily.Isteps, 2);
fnames = fieldnames(traceFamily);
for n = 1 : size(fnames,1)
    param = fnames{n};
    val = getfield(traceFamily, fnames{n});
    %     disp(['size of ' param ' = ' num2str(size(val))])
    if max(size(val)) == 1 % scalars
        try
            disp([prefix fnames{n} ' = ' num2str(getfield(traceFamily, fnames{n})) ';'])
        catch
        end
    elseif max(size(val)) <= numSteps % vectors
        try
            disp([prefix fnames{n} ' = [' num2str(getfield(traceFamily, fnames{n})) '];'])
        catch
        end
    end
end
